clc;clear;clearvars;
% 随机生成10个数据
num_initial = 10;
num_vari = 100;
% 搜索区间
upper_bound = 2.048;
lower_bound = -2.048;
% 迭代时计算y值总的个数
eval_num = 20000;
% K表示划分子空间的个数,sub_num表示每个子空间的维度;k表示领域大小。
S = [5, 10, 25, 50];
K = S(randi(4));
k = 2;
sub_num = num_vari / K;
% 算法的迭代次数,平均子空间的个数为90 / 4 = 22。
K1 = 20;
iter = eval_num / K1;
w = 1;
% 随机生成20个数据,并获得其评估值
sample_x = lhsdesign(num_initial, num_vari).*(upper_bound - lower_bound) + lower_bound.*ones(num_initial, num_vari);
sample_y = Rosenbrock(sample_x);
Fmin = zeros(iter, 1);
aver_Fmin = zeros(iter, 1);

for n = 1 : 100
    n3 = 1;
    % 初始化一些参数
    pbestx = sample_x;
    pbesty = sample_y;
    [fmin, gbest] = min(pbesty);
    global_best_x = pbestx(gbest, :);
    fprintf("n: %.4f\n", n);
    lbest = ones(num_initial, num_vari);
    x = zeros(num_initial, num_vari);
    for i = 1 : iter
        index = randperm(num_vari);
        fming = fmin;
        for i1 = 1 : K
            index1 = 1;
            % 随机分组
            ind = index((1 + (i1 - 1) * sub_num) : i1 * sub_num);
            r = rand(num_initial, sub_num);
            % pso更新下一步的位置,这里可以设置一下超过搜索范围的就设置为边界,服从自由度为1的t分布就是柯西分布
            if r > 0.5
                x(:, ind) = lbest(:, ind) + normrnd(0,1) .* abs(pbestx(:, ind) - lbest(:, ind));
            else
                x(:, ind) = pbestx(:, ind) + trnd(1) .* abs(pbestx(:, ind) - lbest(:, ind));
            end
            x1 = x(:, ind);
            x1(x1 > upper_bound) = upper_bound;
            x1(x1 < lower_bound) = lower_bound;
            x(:, ind) = x1;

            x2 = repmat(global_best_x, num_initial, 1);
            x2(:, ind) = x1;
            y = Rosenbrock(x2);

            % 更新每个单独个体最佳位置
            pbestx1 = pbestx(:, ind);
            pbestx1(y < pbesty, :) = x1(y < pbesty, :);
            pbestx(:, ind) = pbestx1;
            pbesty(y < pbesty, :) = y(y < pbesty, :);
            
            % 更新每个个体的lbest
            n1 = ceil(k / 2);
            n2 = floor(k / 2);
            % 处理开头部分
            for i2 = 1 : n1
                [~, ind1] = min(pbesty(1 : (1 + k), 1));
                lbest(index1, ind) = pbestx(ind1, ind);
                index1 = index1 + 1;
            end
            % 处理中间部分
            for i3 = 1 : (num_initial - k)
                [~, ind2] = min(pbesty(i2 : (i2 + k), 1));
                lbest(index1, ind) = pbestx(i2 + ind2 - 1, ind);
                index1 = index1 + 1;
            end
            % 处理结尾部分
            for i4 = 1 : n2
                [~, ind3] = min(pbesty((num_initial - k) : num_initial, 1));
                lbest(index1, ind) = pbestx(num_initial - k + ind3 - 1, ind);
                index1 = index1 + 1;
            end

            % 更新所有个体最佳位置
            [fmin, gbest] = min(pbesty);
            global_best_x = pbestx(gbest, :);
            %fprintf("iter %d fmin: %.4f\n", i, fmin);
        end
        if fming == fmin
            K = S(randi(4));
            sub_num = num_vari / K;
        end
        Fmin(n3, 1) = fmin;
        n3 = n3 + 1;
    end
    aver_Fmin = aver_Fmin + Fmin;
end
aver_Fmin = aver_Fmin ./ 100;
plot(aver_Fmin);